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Mating preferences of many biological species are not constant but season-dependent. Within the 
framework of evolutionary game theory this can be modeled with two finite opposite-sex populations 
playing against each other following the rules that are periodically changing. By combining Floquet 
theory and the concept of quasi-stationary distributions, we reveal existence of metastable time- 
periodic states in the evolution of finite game-driven populations. The evolutionary Floquet states 
correspond to time-periodic probability flows in the strategy space which cannot be resolved within 
the mean-field framework. The lifetime of metastable Floquet states increases with the size N of 
populations so that they become attractors in the limit N ^ oo. 
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Introduction. The evolutionary dynamics of an animal 
group is tied to the reproductive activity of its mem¬ 
bers, a complex process which involves courtship ritu¬ 
als and sharing of parental care [T]. Within the game 
theory framework, the sex conflict over parental invest¬ 
ment was formalized by Dawkins in his famous “Battle 
of Sexes” (BoS) [2], illustrated in Fig. 1. In this game 
two opposite-sex members of the group play against each 
other. Each player can use two behavioral strategies. 
Entries in the payoff matrix, bss'^ quantify the reward 
received by a female which used a strategy s G {1,2} af¬ 
ter she has played against a male which used a strategy 
s' G {1,2}. Entries Og's define the reward of the male. 

A number of observations have shown that mating 
strategies and preferences of many species are not con¬ 
stant in time but season-dependent [3]. Eor example, 
courtship srituals of the males of Carolina anole lizards 
{Anolis carolinensis), as well as mate selection criteria of 
the females of the species, are periodically changing dur¬ 
ing the year [4]. Even the amount of different types of 
muscle fibers that control the vibrations of a red throat 
fan (dewlap) - which males employ during the courtship - 
is a season dependent characteristic [5] . Currently, there 
is no agreement between the ecologists on the role this 
seasonal plasticity plays in determining the evolution di¬ 
rection of the species [6] . 

We address this problem within the BoS framework 
by allowing the payoffs to periodically vary in time, see 
Eig. 1. Our goal is to investigate how these modulations 
influence the game-driven evolutionary dynamics. Here, 
we first apply the concept of quasi-stationary distribu¬ 
tions in absorbing Markov chains [7] to a stochastic evo¬ 
lutionary dynamics of finite populations and define the 
notion of evolutionary metastable states. Then, by em¬ 
ploying the Eloquet theory we generalize the notion 
of metastable states HDHH] to periodically modulated 



Faithful 


Philanderer 


FIG. 1: (color online) “Battle of Sexes” with seasonal varia¬ 
tions. A female of Carolina anole lizards can be either eoy and 
prefer an arduous courtship, to be sure that a mate is ready 
to contribute to a parental care, or fast, and thus not being 
much concerned about parental care of offspring. A male can 
be either faithful and ready to assure the female partner, by 
performing a long courtship, that he is a faithful potential 
husband, or philanderer and prefer to shorten the courtship 
stage. Depending on the strategies, s {s'), played by the fe¬ 
male (male), the female (male) gets payoff bss' {us's)- Both 
females and males are season-constrained in their strategies 
and preferences, which is modeled via time-periodic modula¬ 
tions of the payoffs. 


game-driven evolutionary dynamics. We show that in 
big but finite populations, the metastable Floquet states 
survive over extremely long (as compared to the period of 
modulations) timescales. We argue that, in the limit of 
infinitely big populations, these states become attractors 
while still evading the mean-field description. 

Model. Finite size of animal populations favors a 





2 


stochastic approach to evolutionary dynamics. Although 
the convergence to the deterministic mean-field dynam¬ 
ics is typically guaranteed in the limit N ^ oc la m, 
the stochastic dynamics of large but finite populations 
can still be very different from the mean-field picture 
[ailBHIH]- Here we adapt the game-oriented version of 
the Moran process m, introduced in Ref. [5] and gen¬ 
eralized to two-player games in Ref. [6]. 

Players A (males) and B (females) form two popula¬ 
tions, each one of a fixed size N and with two avail¬ 
able strategies, s = {1,2}, see Fig. 1. Game payoffs are 
time-periodic functions, Css'{t) = Css'{t + T), c = {a, 6}, 
and can be represented as sums of stationary and zero- 
mean time-periodic components, Css'{t) = Cgs' + Css'{t), 
{css'{t))T = 0. The time starting from t = 0 is incre¬ 
mented by At = T/M after each round. After M rounds 
the payoffs return to their initial values. The state of 
the populations after the m-th round is fully specified by 
the number of players playing the first strategy, i (males) 
and j (females), 0 < i^j < N. A detailed description of 
the corresponding stochastic process is given in Supple¬ 
mental Material. It can also be shown that, in the limit 
^ oc [6], the dynamics of the variables x = i/N and 
y = j/N is defined by the adjusted replicator equations 
[3]; see Refs. [Hin]. 

For a finite N, the state of the system can be expressed 
as a A" X A" matrix p with elements p{i^j)^ which are 
the probabilities to find two populations in the states i 
and jf, respectively. Round-to-round dynamics can be 
evaluated by multiplying the state p with the transi¬ 
tion fourth-order tensor S, with elements 5'(i, j, i', j') [3]. 
By using the bijection k = {N — l)j -h i, we can un¬ 
fold the probability matrix p{i,j) into the vector p(/c), 
/c = 0,..., A^, and the tensor S{i, j') into the matrix 
S{k^l). This reduces the problem to a Markov chain [22] . 
^m+i _ g where m is the number of the round to 
be played. The four states {i = {0, A},j = {0, A}) are 
absorbing states because the transition rates leading out 
of them equal zero [3]. The absorbing states are attrac¬ 
tors of the dynamics for any finite A, and the finite-size 
fluctuations will eventually drive a population to one of 
them [iHl [23] . This would imply a fixation^ so that only 
one strategy survived in each of now monomorphic pop¬ 
ulations [iHH]. 

We are interested in the dynamics before the fixa¬ 
tion, so we merge the four states into a single absorb¬ 
ing state by summing the corresponding incoming rates. 
The boundary states, (i = {0,A},jf G {I,-- - ,A — 1}) 
and (i G {1, • • • , A — 1}, j = {0, A}), can also be merged 
into this absorbing super-state: Once the population gets 
to the boundary, it will only move towards one of the 
two nearest absorbing states. By labeling the absorb¬ 
ing super-state with index /c = 0, we end up with a 


(L + 1) X (T + 1) matrix 


where L = (A — 1)^, is a vector of the incoming 
transition probabilities of the absorbing super-state, 0 is 
a L X 1 zero vector, and is a L x L reduced transition 
matrix. 

With Eq. we arrive at the setup used by Darroch 
and Seneta to formulate the concept of quasi-stationary 
distributions [?]• There is the normalized right eigenvec¬ 
tor of the reduced transition matrix with the maxi¬ 
mum eigenvalue A [24]. By using the inverse bijection, we 
can transform this vector into a two-dimensional proba¬ 
bility density function (pdf), i.e., a state, d, with maxi¬ 
mal mean absorption time. This state is the most resis¬ 
tant to the wash-out by the finite-size fluctuations and it 
remains near invariant, up to a uniform rescaling, under 
the action of the tensor S. This is the metastable state 
of the evolutionary process. 

Stationary case. As an example, we consider a game 
with payoffs an, a22, ^12 and 621 equal 1, and payoffs 
— 1 for the rest of strategies [25]. Figure 2 presents the 
numerically obtained metastable states of the game. We 
use two methods, the direct diagonalization of the re¬ 
duced transition matrix, which is stationary in this case, 
= Q, and preconditioned stochastic sampling [3]. 
For A = 200 we find an agreement between the results 
of the two methods. The means of the metastable state, 

N-l . N-1 . 

i,j=l i,j=l 

coincide with the Nash equilibrium [26] for any A. 
However, the actual dynamics is determined by the 
metastable limit cycle encircling the equilibrium (this 
could be seen by performing short-run stochastic sim¬ 
ulations); see Fig. 2. Within the Langevin-oriented 
approach to the dynamics of finite populations 
the appearance of the metastable limit cycle can be in¬ 
terpreted as a stochastic Hopf bifurcation [28] (see also 
Ref. [29] for another interpretation). In the limit A ^ 00 
the cycle collapses to the Nash equilibrium. Note, how¬ 
ever, that the convergence to this limit is slow, as indi¬ 
cated by the width of the pdf for A = 400. 

Case of modulated payoffs. By adding time- 
modulations to the model, we find that the mean-field 
dynamics does not exhibit substantial changes. For 
the choice e{t) = aii(t) = 622 (^) = f cos{ujt) with 
uj = 27t/T (all other payoffs held stationary) we ob¬ 
served a period-one limit cycle localized near the Nash 
equilibrium of the stationary case, see Figs. 3(a,b). 
It collapses to a set of adiabatic Nash equlibria, 
|a;jv_E(e) = |Ef, 2/JV_E(e) = in the limit w 0. 
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FIG. 2: (color online) Metastable states of the stationary BoS game. In the mean-filed limit A/" = oo, a trajectory spirals towards 
a fixed point (|, |), the Nash equilibrium of the game. For the finite N, metastable states are specified by their quasi-stationary 
probability density functions (pdf’s) (3d plots). For N = 200 the pdf combines the results of the direct diagonalization of the 
39 601 X 39 601 matrix Q (left half of the pdf, this procedure was also used to obtain the function for N = 100) and of the 
preconditioned stochastic sampling (right part of the pdf, this procedure was also used to obtain the function for N = 400) [3]. 
The baseline fitness w = 0.3 (other parameters are given in the text). 


The dynamics of a finite N population is different. The 
stochastic evolution of a trajectory in the (i, j)-space, ini¬ 
tiated away from the absorbing boundary, can be divided 
into two stages. At first the trajectory relaxes towards 
a metastable state. The timescale of this process is de¬ 
fined by the mixing time [30], which in this case 

has to be calculated now for the g^/asTstationary state. 
Then the trajectory wiggles around the metastable state 
until the fluctuations drive it to the absorbing boundary. 
Following the random-walk approximation, the mean ab¬ 
sorption time tabs (A^) 5 called “mean fixation time” [ 2 l[TT] 
in the evolutionary context, seemingly should also scale 
as N. However, this estimate neglects the presence 
of the inner attractive manifold and the fact that the 
noise strength decreases upon approaching the absorb¬ 
ing boundary. In fact, the absorption time scales super- 
linearly with N 1ST]. The lifetime of the metastable state 
is restricted to the time interval [tmix(A^), ^abs(A^)]5 whose 
length scales as 4bs(V)[l - W(V)/tabs(V)] ~ tabs(V). 

For u) = 0.1 |32], the stochastic simulations reveal 
a metastable state which is distinctively different from 
the limit cycle produced by the mean-field equations, 
see Fig. 3b. There is a conflict between the evolution 
of means, described by the adjusted replicator equa¬ 
tions, and the results of the stochastic dynamics. The 
conflict can be resolved with the concept of the quasi¬ 
stationary distribution. Namely, the transition matrices, 
Eq. Q, are round-specific now and form a set {S"^}, 
m = 1, • • • , M. The propagator over the time interval 
[0,t], 0 < t < T, is the product U(t) = Ylm^=i with 
Mf = t/At. All the propagators, including the period- 
one propagator U(T), have the same structure as the 
super-matrix in Eq. 0- We define the metastable state 
d(T) as the the quasi-stationary distribution of U(T). 
It is also a Floquet state |9] of the reduced propagator 
U”(T), which can be obtained by replacing the tran¬ 
sition matrices with the matrices or by sim¬ 

ply cutting out the first line and column from the ma¬ 


trix U(T). The Eloquet state is a time-periodic state, 
d(t + T) = d(t), which changes during one period of 
modulations, see Eig. 4. The metastable state d{t) at 
any instant of time t, 0 < t < T, can be obtained by act¬ 
ing on the state d(0) with the reduced propagator U^(t). 

The evolution of the means of the pdf d(t) (see Eig. 4a), 
(T(t),^(t)), is close to the period-one limit cycle, see blue 
dots on Eig. 3a. However, the Eloquet state consists of 
two peaks produced by the noised period-two limit cycle 
(compare also the positions of the stroboscopic points 
in Eig. 3b with the pdf for t = 0 in Eig. 4a). The 
peak contributions balance each other thus reducing the 
dynamics of the means to the vicinity of the the point 
(^, ^). The lifetime of the state d(t) can be estimated 
with the largest eigenvalue At, 0 < At < 1, of the ma¬ 
trix U’^(T). To compare it with lifetimes of stationary 
metastable states, we introduce the mean single-round 
exponent. At = A^^^ and define the mean lifetime as 
^life = 1/(1 ~ ^t) [3|. Aside of the slow decay trend, we 
found the effect of modulations not being strong. This 
is in stark contrast to the structure of the metastable 
states. Namely, while in the stationary limit the pdf d 
is localized near the Nash equilibrium, at the maximal 
distance from the absorbing boundaries, the metastable 
Eloquet state is localized near the absorbing boundary, 
see Eig. 4. We also detect the increase of the boundary 
localization with the increase of the population size be¬ 
yond N = 200. This suggests that, in the limit A ^ oo, 
the dynamics of the system is governed by a period-two 
limit cycle localized near the absorbing boundary. The 
boundary localization of the metastable attractor can be 
interpreted as the presence of small fractions of mutants 
in, i.e. the players that are using strategies different 
from that used by the majority of populations. The evo¬ 
lutionary dynamics of the mutant fractions looks like a 
repeating sequence of population bottlenecks [U [33] yet 
this only weakly affects fraction lifetimes [34] even in the 
case of finite N. 
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FIG. 3: (color online) Evolutionary dynamics governed by the 
BoS game with modulated payoffs, (a) Period-one limit cy¬ 
cles of the mean-field dynamics for ct; = 0.1 (blue dash-dotted 
line) and uj = 0.01 (red solid line) are localized near the Nash 
equilibrium of the stationary game, (|, |) (arrows indicate 
the direction of motion). In the limit cj ^ 0, the mean-field 
attractor shrinks to the set of adiabatic Nash equlibria (black 
dashed line). Mean position (•) of a finite-iV metastable Flo- 
quet state, Eq. 0, moves along the limit cycle 

localized near the point (|, fj (the means are plotted at the 
instants tn = nT/5, n = 0, ..,4); (b) A stochastic trajectory 
(grey line) reveals the existence of a period-two limit cycle 
[the period doubling can be resolved with stroboscopic points, 
plotted at the instants 2nT (A) and (2n+ 1)T (<C>)]. The tra¬ 
jectory is initiated at the point marked with the open blue 
square and ends up at the absorbing state (red cross at the 
upper left corner). The trajectory of the mean of the finite- 
N metastable Floquet state (•) is distinctively different from 
the stochastic trajectory [note the change of scale as com¬ 
pared to panel (a)]. The parameters are / = 0.5, N = 200, 
and M — T!l\t — lOiV (corresponds to the driving frequency 
cc = 0.1 in the mean-field limit) [32]. Other parameters as in 
Fig. 2. 


Conclusions. We presented a concept of metastable 
Floquet states in game-driven populations when mate 
selection preferences are periodically changing in time. 
Here we combined the Floquet formalism with the con¬ 
cept of quasi-stationary distributions to reveal the exis¬ 
tence of complex, liquid-like nonequlibrium dynamics in 
the strategy space which cannot be resolved within the 
mean-field framework. Metastable Floquet states are not 
restricted to the field of ecology studies but can emerge in 
different periodically modulated systems with stochastic 
event-driven dynamics. They may, for example, underlay 
a gene expression in a single cell, which is modulated by 
a circadian rhythm [38] and can provide new interpre¬ 
tations of the Bose-Einstein condensation in ac-driven 
atomic ensembles [3S1I10]. 
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Supplemental Material 


SEASONAL VARIATIONS IN MATE 
PREFERENCES 

Stable time variations were found in the female fly¬ 
catcher preferences or male forehead patch size that 
resulted in late-breeding females preferring males with 
larger patches [Tj. It was explained by the fact that in 
the beginning of the breeding season, large-patched males 
allocate more resources to courting than to parental care 
but change their habits to the opposite late in the season. 
Seasonal variations were also found in fiddle crabs (female 
preference to male claw size) [2], two-spotted goby (fe¬ 
male preference to overall male size) [3], and sailfin mol¬ 


lies (male preferences for two different kind of females) 

a- 


MORAN PROCESS 


Players A (males) and B (females) form two popula¬ 
tions, each one of a fixed size N and with two avaiiabie 
strategies, s = {1,2}. Payoffs are specified by four func¬ 
tions, {ass'{t)} and {65/5(t)}, 5,5' = {1,2}. The average 
payoff of the piayers using strategy s is 

EU^t)=asi{t)j^+as2{t)SE^, (3) 

Trf{i,t) = bsi{t)^ + bs2{t)SL^_ (4) 


Payoffs determine the probabiiities for a piayer to be cho¬ 
sen for reproduction, e.g. for the maie popuiation. 


^ ’ N 1 — re + re7r^(i, jf, t) ’ 


(5) 


where S- {N — 'i)7r^(jf, t)]/V is the 

average payoff of the maies. The baseiine fitness w G 
[0,1] is a tunabie baseiine fitness parameter determining 
how the piayer’s chance to be chosen for reproduction 
is reiated to piayer’s performance [Si]. When re = 0, 
the probabiiity to be chosen for reproduction does not 
depend on piayer’s performance and is uniform across 
the popuiation. 

After the choice has been made, another member of the 
popuiation is chosen compieteiy randomiy and repiaced 
with an offspring of the piayer chosen for reproduction, 
i.e. with a piayer using the same strategy as its par¬ 
ent [7]. This update mechanism is acting simuitaneousiy 
in both popuiations, A and 5, such that a mating pair 
produces two offspring, a maie and a femaie, on every 
round. Therefore, the size of the popuiations N remains 
constant. 

A single round can be considered as a one-step Markov 
process, with transition rates, e.g. for population A, from 
a state i to states i + 1 and i — 1, are given by mun] 


TA(.hj,t) 




1 — w S-WTTi (t) i N — i 
1 — w W7t^ N N 
1 — w S- WTT^it) N — i i 
1 — W WTT^ N N 


( 6 ) 


TRANSITION TENSOR 

Here we describe the transition fourth-order tensor 
jA'^ j') in terms of the rates [T^’~(i, j, t) and 
Tb'~ ih jA)] for populations A and B given by Eq. (4) 
in the main text. The stochastic Moran process can be 
expressed as a Markov chain m 
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= [l-TX{i,j:m/S.t) -T^(i,j,mM)\ [l-T^{i,j,m/S.t)-Tg{i,j,m/S.t)\p^(i,j) 

+ l,TOAt) [l - + l,mAi) - TX{i,j + p^{i,j + 1) 

+^"b(*A' - l,mAt) [1 -T^(z,j - l,mAt) - l,mAt)]p'^{i,j - 1) 

+T^{i + l,j,mAt) [1 -Tg{i + l,j,mAt) -T^{i + l,j,mAt)] p^{i + l,j) 

+TX (* - 1, j, mAt) [1 l,j, mAt) l,j, mAt)] p'^ii - 1, j) 

+T^ (i + 1, jf + 1, m/\t)T^ (i + 1, j + 1, m/\t)p^ [i H- 1, j + 1) 

+Tpi - 1, j + l,mAi)T3 (z - 1, j + l,mAt)p"^{i -l,j + 1) 

+TX{i + l,j - l,mAt)Tpi + l,j - l,mAt)p"^{i + l,j - 1) 

+Tpi -l,j - l,mAt)Tpi -l,j - l,mAt)p^{i -l,j - 1). (7) 


The above equation can be recast into where the fourth-order tensor S^{i, , j') is given by, 


= [l-T+{i',j',mAt)-TX{i',f,mAt)] [l - T+{i',j',mAt) - T^ii',j',mAt)] 6i,,i6pj 
+^B(*^i^”^A^) [l -TXii',j',mAt) -T^{i',j',mAt)] 5i>^i5pj+i 
+TB{i',j',mAt) [l - TX{i',j\mAt) - T^{i',j',mAt)] 

+TX{i',j','rnAt) [l - Tg {i', j', mAt) - {i', j', mAt)] Si>^i+i5pj 

+^4 (i',/,TOAi) [l - TB{i',j\mAt) - T^{i',j',mAt)] 

+T^mAt)TB wAi) V,i+i 

+Ta mAt)TB mAt)6i,^i_i 5pj+i 

+TX mAt)T^mAt)Si>^i+i Sj> j-i 

+Ta mAt)T^ (z',/, mAt)Si>^i-i Spj-i. (9) 


Above i = 0, • • • , A^, j = 0, • • • , A^, i' = 0, • • • , A^, and 
/ = 0, • • • ,N. Using the bijection k = {N — l)j + i and 
I = [N — l)j' + i', we obtain the required matrix form, 
see Eq. (7) in the main text. 

ADJUSTED REPLICATOR EQUATIONS 

In the continuous limit A ^ oc, the dynamics of the 
variables x = i/N and y = j/N is defined by the adjusted 
replicator equations mill], 

( 10 ) 

^ = + (11) 

where = C12-C22, = C11+C22-C12-C21, T = 

and C = {A^B}. Tt^(x,y,t) [tt^{ x^y,t)\ is the averaged 
(over the population) payoff of the males [females]. 


THE LIFETIME OF A METASTABLE STATE 

The lifetime of the state d(t) can be estimated with 
the largest eigenvalue At, 0 < At < 1, of the ma¬ 
trix U’^(T). To compare it with lifetimes of stationary 
metastable states, we introduce the mean single-round 
exponent. At = A^^^ and define the mean lifetime as 
hife = 1/(1 — At). Figure 1 shows the dependence of tufe 
on the strength of modulations. 

SIMULATIONS 

The preconditioned stochastic sampling was performed 
by launching trajectories from random initial points, uni¬ 
formly distributed on the A — 1 x A — 1 grid and then 
sampling the pdf with only those trajectories which re¬ 
mained unabsorbed after 10 • A^ rounds. 

For A = 200 the diagonalization of the 39 601 x 39 601 
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FIG. 5: (color online) The lifetime tufe as a function of the 
modulation strength /, for the population size N = 50 (O)^ 
100 (□), and 200 (A). Other parameters are as in Fig. 2 in 
the main text. 

matrix Qt was performed on the cluster of the MPIPKS 
(Dresden) and Leibniz-Rechenzentrum (Miinchen). The 
stochastic sampling was performed on a GPU cluster con¬ 
sisting of twelve TESLA K20XM cards. That allowed us 
to obtain 5 • 10^ realizations for each set of parameters. 
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